Orthogonal Polynomial Representation of Imaginary-Time Green's Functions 
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We study the expansion of single-particle and two-particle imaginary-time Matsubara Green's 
functions of quantum impurity models in the basis of Legendre orthogonal polynomials. We discuss 
various applications within the dynamical mean-field theory (DMFT) framework. The method pro- 
vides a more compact representation of the Green's functions than standard Matsubara frequencies 
and therefore significantly reduces the memory-storage size of these quantities. Moreover, it can be 
used as an efficient noise filter for various physical quantities within the continuous-time quantum 
Monte Carlo impurity solvers recently developed for DMFT and its extensions. In particular, we 
show how to use it for the computation of energies in the context of realistic DMFT calculations in 
combination with the local density approximation to the density functional theory (LDA+DMFT) 
and for the calculation of lattice susceptibilities from the local irreducible vertex function. 

PACS numbers: 71.27.+a, 71.10.Fd 



In recent years, significant progress has been made in 
the study of strongly-correlated fermionic quantum sys- 
tems with the development of methods combining sys- 
tematic analytical approximations and modern numerical 
algorithms. The Dynamical Mean-Field Theory (DMFT) 
(for a review see Ref. 1) and its various extensions. 2-6 
serve as successful examples for this theoretical advance. 
On the technical side, important progress was made in 
the solution of quantum impurity problems, i.e. local 
quantum systems coupled to a bath (self-consistently 
determined in the DMFT formalism). In particular, a 
new generation of continuous-time quantum Monte Carlo 
(CTQMC) impurity solvers 7-10 has emerged that provide 
unprecedented efficiency and accuracy (for a recent re- 
view, see Ref. 11). 

In practice, several important technical issues still re- 
main. Firstly, while the original DMFT formalism is 
expressed in terms of single-particle quantities (Green's 
function and self-energy), two-particle quantities play a 
central role in the formulation of some DMFT extensions 
(e.g. dual-fermions 4,12-14 , DrA 3 ) and in susceptibility 
and transport computations in DMFT itself. They typi- 
cally depend on three independent times or frequencies, 
and spatial indices. Therefore, they are quite large ob- 
jects that are hard to store, manipulate and analyze, even 
with modern computing capabilities. Developing more 
compact representations of these objects and using them 
to solve, e.g., the Bethe-Salpeter equations is therefore 
an important challenge. 

A natural route is to use an orthogonal polynomial 
representation of the imaginary-time dependence of these 
objects. While the application of orthogonal polynomials 
has had productive use in other approaches to correlated 
electrons, 15 ' 16 in this paper we show how to use Legendre 
polynomials to represent various imaginary-time Green's 
functions in a more compact way and show their useful- 
ness in some concrete calculations. 

A second aspect is that modern CTQMC impurity 
solvers still have limitations. One well-known problem 



is the high-frequency noise observed in the Green's func- 
tion and the self-energy (see e.g. Fig. 6 of Ref. 17). 
Even though this is in general of little concern for the 
DMFT self-consistency itself, it can become problematic 
when computing the energy, since the precision depends 
crucially on the high-frequency expansion coefficients of 
the Green's function and self-energy. An important field 
of application involves realistic models of strongly cor- 
related materials through the combination with the lo- 
cal density approximation (LDA+DMFT). 18 In this pa- 
per, we show that physical quantities such as the Green's 
function, kinetic energy, and even the coefficients of the 
high-frequency expansion of the Green's function can be 
measured directly in the Legendre representation within 
CTQMC and that the basis truncation acts as a very ef- 
ficient noise filter: the statistical noise is mostly carried 
by high-order Legendre coefficients, while the physical 
properties are determined by the low-order coefficients. 

This paper is organized as follows: Section I is devoted 
to single-particle Green's functions. More precisely, in 
Sec. LA, we introduce the Legendre representation of 
the single-particle Green's function and how it appears 
in the CTQMC context; we then illustrate the method on 
the imaginary-time (LB) and imaginary-frequency (I.C) 
Green's function of a standard DMFT computation; in 
I.D, we discuss the use of the Legendre representation to 
compute the energy in a realistic computation for SrV03. 
Section II is devoted to two-particle Green's functions: 
We first present the expansion in Sec. II.A and illustrate 
it on an explicit DMFT computation of the antiferromag- 
netic susceptibility in Sec. II. B, followed by the example 
of a calculation of the dynamical wave-vector resolved 
magnetic susceptibility. Additional information can be 
found in the appendixes. Appendix A gives some prop- 
erties of the Legendre polynomials relevant for this work. 
Appendix B discusses the rapid decay of the Legendre 
coefficients of the single-particle Green's function. Ap- 
pendix C first derives the accumulation formulas for the 
single-particle and two-particle Green's functions in the 
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hybridization expansion CTQMC (CT-HYB) algorithm 8 
(while these formulas have been given before, 8,17 the 
proof presented here aims to explain their resemblance to 
a Wick's theorem). We then give the explicit formulas in 
the Legendre basis. For completeness, we provide an ac- 
cumulation formula for the continuous- time interaction 
expansion (CT-INT) 7 and auxiliary field (CT-AUX) 10 al- 
gorithms in Appendix D. Finally, in Appendix E, we 
derive the expression for the matrix that relates the co- 
efficients of the Green's function in the Legendre repre- 
sentation to its Matsubara frequency representation. 

I. SINGLE-PARTICLE GREEN'S FUNCTION 

A. Legendre representation 

We consider the single-particle imaginary-time Green's 
function G(r) defined on the interval [0, /3], where j3 is the 
inverse temperature. Expanding G(t) in terms of Legen- 
dre polynomials Pi{x) defined on the interval [—1,1], we 
have 



i>o P 

G l =V2l + T <1tPi(x{t))G{t). (2) 
Jo 

where x(r) = 2r//3 — 1 and Gi denote the coefficients of 
G(t) in the Legendre basis. The most important prop- 
erties of the Legendre polynomials are summarized in 
Appendix A. 

We note that a priori different orthogonal polynomial 
bases (e.g., Chebyshev instead of Legendre polynomials) 
may be used, and many of the conclusions in this pa- 
per would remain valid. The advantage of the Legendre 
polynomials is that the transformation between the Leg- 
endre representation and the Matsubara representation 
can be written in terms of a unitary matrix, since Legen- 
dre polynomials are orthogonal with respect to a scalar 
product that does not involve a weight function (see be- 
low and Appendix E). In this paper, therefore we restrict 
our discussion to the Legendre polynomials. 

On general grounds, one can expect the Legendre rep- 
resentation of G(t) to be much more compact than the 
standard Matsubara representation: in order to perform 
a Fourier series expansion in terms of Matsubara frequen- 
cies, G(t) has to be anti-periodized for all r € R, while 
the full information is already contained in the interval 
[0, /?]. As a result, the Green's function contains disconti- 
nuities in t that result in a slow decay at large frequencies 
(typically ~ \jv n ). On the other hand, expanding G(t), 
which is a smooth function of r on the interval [0, /3], in 
terms of Legendre polynomials yields coefficients Gi that 
decay faster than the inverse of any power of I (as shown 
in Appendix B). As a result, the information about a 
Green's function can be saved in a very small storage 
volume. As we will show in Sec. II, this is particularly 
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FIG. 1. (Color online) Legendre coefficients Gi of the Green's 
function of the half-filled Hubbard model on the Bethe lattice 
within DMFT. Error bars are not shown on the logarithmic 
plot. They are of the order of 10 -4 . 

relevant when dealing with more complex objects such 
as the two-particle Green's function, which depends on 
three frequencies. 

CTQMC algorithms usually measure the Green's func- 
tion G(t) in one of the two following ways: (i) using a 
very fine grid for the interval [0,/3] or (ii) measuring the 
Fourier transform of the Green's function on a finite set 
of Matsubara frequencies. 7,19 We show in Appendix C 
cxplicitcly for the CT-HYB 8,11 algorithm, that one can 
also directly measure the coefficients Gi during the Monte 
Carlo process (we expect our conclusions to hold for any 
continuous-time Monte Carlo algorithm). 

As an illustration, we will focus on the Green's function 
obtained by DMFT for the Hubbard model at half-filling 
described by the Hamiltonian 

H = -t 4* c j<? + Uj2n it n 4 , (3) 

(ij)<7 i 

where cfj creates (annihilates) an electron with spin a 
on the site i of a Bcthc lattice 1,20 and (ij) on the sum 
denotes nearest neighbors. In the following, quantities 
will be expressed in units of the hopping t and we set 
the on-site Coulomb repulsion to U/t = 4 and use the 
temperature T/t = 1/45. We solve the DMFT equations 
using the TRIQS 21 toolkit and its implementation of the 
CT-HYB 8,11 algorithm. In Fig. 1, we show the coeffi- 
cients Gi that we obtain. Note that coefficients for I odd 
must be zero due to particle-hole symmetry. Indeed, the 
coefficients in our data for odd Vs all take on very small 
value, compatible with a vanishing value within their er- 
ror bars. The even I coefficients instead show a very fast 
decay, as discussed above. For I > 30, all coefficients 
eventually take values of the order of the statistical error 
bar. 

Let us now discuss the specific issue of the statistical 
Monte Carlo noise. We observe that the high-order Leg- 
endre coefficients have a larger relative noise than small 
I coefficients. On general grounds, we expect the co- 
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FIG. 2. (Color online) Imaginary-time Green's function G(r) 
at four different values of r as a function of i max - 



efficients of the exact Green's function to continue to 
decrease faster than any power of l/l to zero (cf. Ap- 
pendix B). Hence, physical quantities computed from 
G(r) arc likely to have a very weak dependence on the 
G[ for large I. A good approximation then is to trun- 
cate the expansion in Legendre polynomials at an order 
? max and set Gi = for / > Z max . The choice for Z max 
has to be such that the quantity of interest is accurately 
represented. On the other hand, if Z max is too large, we 
would start to include coefficients that have increasingly 
large error bars compared to their value and this would 
eventually pollute the calculation. A systematic method 
is therefore to examine the physical quantity as a func- 
tion of the cutoff / max - We expect that it will first reach 
a plateau where it is well converged. The existence of 
a plateau means that the contribution of higher-order 
coefficients is indeed negligible. For larger Z max , the sta- 
tistical noise in the Gi will destabilize this plateau whose 
size will increase with the precision of the CTQMC com- 
putation. The existence of such a plateau provides a 
controlled way to determine the adequate value of Z max . 
Tn the remaining paragraphs of this section, we will il- 
lustrate this phenomenon on different physical quantities 
by studying their dependence on i max . 



FIG. 3. (Color online) Imaginary-time Green's function G(r) 
on the interval [0, /3] measured on a finite 1500-bin mesh (blue 
scattered points) and computed from Z max Legendre coeffi- 
cients (solid lines). Four different choices for i max are shown. 
Inset: zoom on the area around (5/2. 

values of G(r) have not yet converged to their plateau 
(see Fig. 2), the resulting Green's function is smooth but 
not compatible with the scattered direct measurements. 
For l max = 35 and 60, G(t) is smooth and nicely in- 
terpolates the scattered data. Moreover G(t) is virtually 
identical for both values of Z max . This is expected because 
both of these values lie on the plateau. When / max is very 
large, i.e., of the order of the number of imaginary-time 
bins, the noise in G(r) eventually reappears and begins to 
resemble that of the direct measurement. We emphasize 
that all measurements have been performed within the 
same calculation and hence contain identical statistics. 
Hence the information in both measurements is identical 
up to the error committed by truncating the basis. 

It is clear from this analysis, that the truncation of 
the Legendre basis acts as a noise filter. We note that 
no information is lost by the truncation: the high-order 
coefficients correspond to information on very fine details 
of the Green's function, which cannot be resolved within 
a Monte Carlo calculation, as is obvious from the noisy 
G(r). 



B. Imaginary-time Green's function 

It is instructive to analyze the effect of Zmax on the 
reconstructed imaginary-time Green's function G(r) (us- 
ing Eq. (1)). In Fig. 2, wc show the evolution of G(t) at 
t = 0+, r = f3/8, t = p/4, and r = /3/2 with the cutoff. 
It is apparent that these values very rapidly converge as a 
function of Z m ax- We observe a well-defined and extended 
plateau. As the cutoff grows bigger, noise reappears in 
G(t) because of the comparatively large error bars in 
higher-order G;'s. 

In Fig. 3, the Green's function is reconstructed on the 
full interval [0, /?] and compared to a direct measurement 
on a 1500-bin mesh. For Z max = 20, where the individual 



C. Matsubara Green's function and high-frequency 
expansion 

It is common to use the Fourier transform G(i^ n ) of 
G(r) to manipulate Green's functions. This represen- 
tation is, for example, convenient to compute the self- 
energy from Dyson's equation or to compute correlation 
energies. In terms of G;, we can obtain the Matsubara 
Green's function with 

G{\p n ) =VG,^ [ dr e^ T fl(x(r)) 
;>o P Jo 

= J2 T "iGi. (4) 
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FIG. 4. (Color online) Matsubara Green's function obtained 
from measurements made directly on the Matsubara frequen- 
cies (blue scattered points), calculated from an imaginary- 
time measurement (green scattered points) and computed 
from Eq. (4) with Z max = 35 (red solid line). The ana- 
lytically known high-frequency tail is shown for comparison 
(black solid line). Inset: Blowup of the high-frequency region. 



where the unitary transformation T n i is shown in Ap- 
pendix E to be 

^^(-ir^V^TTJ^^-) (5) 



with ji(z) denoting the spherical Bessel functions. Note 
that T n i is independent of /3. 

In Fig. 4, we display the Matsubara Green's function 
as measured directly on the Matsubara axis and as com- 
puted from Eq. (4) with a fixed cutoff Z max . The direct 
measurement of G{\v n ) has been done within the same 
Monte Carlo simulation as the one used to compute the 
Gi discussed above. It is clear from the plot that the trun- 
cation to Z max has filtered the high-frequency noise, and 
that for large \v n the Matsubara Green's function has a 
smooth power-law decay. Let us emphasize here that the 
Matsubara Green's function is obtained in an unbiased 
manner that does not involve any model-guided Fourier 
transform (see also Ref. 22). 

We will now show that the coefficients that control this 
power-law decay can also be accurately computed. Let 
us consider the high-frequency expansion of G(\v n ) 



G(iu n ) = A- 



C2 



<:':! 



(iv n ) 2 O'n) 3 



(6) 



Using the known high-frequency expansion of T„; (cf. 
Appendix E), 



,(i) 



.(2) 



f (3) 



T n l — — 



(7) 



one can directly relate the c v and the Gi. Indeed, from 
(4), (6), and (7), it follows that 



(8) 



FIG. 5. (Color online) Convergence of the moments ci, cz 
and C5 as a function of Z max . Only points corresponding to an 
even cutoff are shown because odd terms in the sum vanish. 
The analytically known results for c\ and C3 are indicated 
by dashed lines. Even moments are zero due to particle-hole 
symmetry. 



The general expression of the coefficients tj is shown 
in (E8). For the first three moments, we have the follow- 
ing expressions 



ci 



C2 = 



C3 = 



;>o, 



;>0, odd 



(9a) 
(9b) 



E ^T~ G 



Z>0, even 



Since ~ ;2p-3/2 j with the fast decay of the G[ dis _ 
cussed above, we can expect a stable convergence of the 
c p as a function of Z max - Note, however, that when p in- 
creases, the coefficients grow, so we expect to need more 
and more Legendre coefficients to compute the series in 
practice. 

The convergence of the moments is illustrated in Fig. 5. 
For the model we consider, the first moments are explic- 
itly given by 

C\ = 1 C2 = 

C3 = 5 C4 = 0. 

We see that c\ and C3 smoothly converge to a plateau. 
For the higher moment C5, a larger number of Legen- 
dre coefficients is required. A plateau is reached but is 
(depending on the accuracy of the data accumulated in 
the QMC simulation) quickly destabilized when Z max gets 
bigger and noisy Gi are included in the calculation. This 
clearly shows that Z max has to be chosen carefully to get 
sensitive c p . For larger cutoff the error in the moments 
grows rapidly. This shows that a large error on the high- 
frequency moments is committed when measuring in a 
basis in which it is not possible to filter the noise, i.e., 
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the conventional imaginary-time or Matsubara represen- 
tation. 

Note that it is easy to incorporate a priori information 
on the moments c p . For example, in the model we con- 
sider above, we have ci = ({c, c 1 *}) = 1. From (9a), we 
see that this is a linear constraint on the Gi coefficients, 
which we can therefore enforce by projecting the Legen- 
dre coefficients onto the (Z max )-dimcnsional hyperplanc 
defined by the constraint (9a). A correction to impose, 
e.g., a particular c\ is straightforwardly found to be 

Gi^G l+ (fa E^G.) -f^y-. (10) 
\ i'=o / l^i \h I 

This is easily generalized to other constraints. 

D. Energy 

The accurate determination of the high-frequency co- 
efficients is of central importance, since many quantities 
are computed from sums over all Matsubara frequen- 
cies involving G(iu n ). Because G(iv n ) slowly decreases 
as ~ l/(iv n ) to leading order, these sums are usually 
computed from the actual data up to a given Matsubara 
frequency and the remaining frequencies are summed up 
analytically from the knowledge of the c p . Thus, an incor- 
rect determination of the c p leads to significant numerical 
errors. This is a particularly delicate issue when G(w n ) 
is measured directly on the Matsubara axis. In this case 
one usually needs to fit the noisy high-frequency data to 
infer the high-frequency moments. As discussed above, 
such a procedure is not required when using Legendre 
coefficients and the c p can be computed in a controlled 
manner. In the following, we illustrate this point in an 
actual energy calculation. 

Based on an LDA+DMFT calculation for the com- 
pound SrVC>3, 23,24 we compute the kinetic energy E^in = 
{l/N)J2 k . a ( n ka}£ka and the correlation energy E co „ = 
Si Uin^rii]) (N denotes the number of lattice 
sites) resulting from the implementation and parame- 
ters of Rcf. 24. These terms are contributions to the 
LDA+DMFT total energy 25 which depend explicitly on 
the results of the DMFT impurity solver. 

The results are shown in Fig. 6. Here the parameter 
'rnaxj against which these quantities are plotted, repre- 
sents the number of Legendre coefficients used through- 
out the LDA+DMFT self-consistency. It is also the num- 
ber of coefficients used to evaluate (rika) from the lattice 
Green's function Gfe(if n ). Note that S corr has been ac- 
cumulated directly within the CTQMC simulation. 

In agreement with an analysis of the convergence with 
respect to the number of Legendre coefficients l max sim- 
ilar to the ones shown in Figs. 2, 5 for an individual 
DMFT iteration, we find a plateau for both energies at 
( max ~ 40. While the energy can be accurately computed 
within a single DMFT iteration, the error here mainly 
stems from the fluctuations between successive DMFT 
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FIG. 6. (Color online) Kinetic energy _Eki n (full symbols) and 
correlation energy -E CO rr (open symbols) for SrVC>3 as a func- 
tion of imax, computed with the implementation and parame- 
ters of Ref. 24. For clarity the kinetic energy has been shifted 
by 384.86eV. Error bars are computed from 80 converged 
LDA+DMFT iterations. 

iterations. The plateau remains up to the largest values 
of ?max- However, as Z max gets larger, so do the error bars, 
due to the feedback of noise from the largest Legendre 
coefficients. Note that the error bars on the correlation 
energy, computed directly within the CTQMC algorithm, 
are of the same order of magnitude as those on the ki- 
netic energy. The existence of a plateau implies that for 
a well-chosen cutoff / ma x, the energy can be computed in 
a controlled manner. We want to emphasize that such an 
approach is simpler and better controlled than delicate 
fitting procedures of high-frequency tails of the Green's 
function on the Matsubara axis. 



II. TWO-PARTICLE GREEN'S FUNCTION 

A. Legendre representation for two-particle 
Green's functions 

The use of Legendre polynomials proves very useful 
when dealing with two-particle Green's functions. We 
will show that it brings about improvements both from 
the perspective of storage size and convergence as a func- 
tion of the truncation. The object one mainly deals with 
is the generalized susceptibility, 

(n2, T34, T M ) = X a<J ( T l - T 2, 73 ~ T 4 , T\ - T 4 ) = 
(T4(Vl )Ca(T 2 )c\, (t 3 )c^(t 4 )) 

- (TcUnfafaMT&inMTj). (11) 

Let us emphasize that x is a function of three indepen- 
dent time-differences only. With the particular choice 
made above, \ is /3-antiperiodic in t\% and T34, while it 
is /3-periodic in T14. Consequently, its Fourier transform 
x(iv ni \v n i , iw m ) is a function of two fermionic frequencies 
v n = 2(n + 1)tt//3, v n i = 2(n' + l)ir//3, and one bosonic 
frequency uj m = 2m7r//3. 
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We introduce a representation of x( T i2, T 34, T u) in 
terms of the coefficients Xll'i^m) such that 



X(T12,T34,T 14 ) = 2^ 2^ 03 

(,l'>Om£Z 

J Pi(^(n2))P/'(x(r 3 4))e iw '" Tl4 xzr(iwm). 



f 



(-1)' 



112) 



In this mixed basis representation, the t\i and T34 de- 
pendence of x( T i2 7 T 34, T i4) is expanded in terms of Leg- 
endre polynomials, while the T14 dependence is described 
through Fourier modes e 1WmTl4 . The motivation behind 
this choice is that many equations involving generalized 
susceptibilities (like the Bcthe-Salpetcr equation) are di- 
agonal in \uj m . The inverse of (12) reads 



dT 12 dT 34 dT 1 W'2l + ly/21' + 1(-1)'' + 1 



Pi (x(t 12 ))P v (x(r 34 ))e- i ^^x(Ti2, r 34 , r 14 ). 

(13) 



We show in Appendix C how the Legendre expansion 
coefficients of the one- and two-particle Green's function 
(hence of xiv (iw m ) ) can be measured directly within CT- 
HYB. With the above definition, the Fourier transform 
x(iv n , 'w n ', iw m ) is easily found with 



X(iz^n,hv,iw m ) = 2J T n ixii'{iu m )T* n ,. 



(14) 



l,l>0 



T n i was already defined in Eq. (4). Using the additional 
unitarity property of T in Eq. (14) one can in general 
easily rewrite equations involving the Fourier coefficients 
x(iv n , il V, iw m ) in sole terms of the x;z< (iw m ). 

In the DMFT framework, the lattice susceptibility xiatt 
is obtained from 1,2 



Xlatt 



(i^m,q) 



\lo 



Ylo 



-1 

(iWm) 

-1 



V° 
Alatt 



(iw 



(15) 



where the double underline emphasizes that this is to be 
thought of as a matrix equation for the coefficients x ex- 
pressed either in (iu n , \v n i ) in the Fourier representation 
or in (I, I') in the mixed Legendre- Fourier representation. 
The bare susceptibilities are given by 

Xi oc (nv„,iz/„/,ia; m ) = -Gi oc (w n + iu} m )G\ oc (iv n )5 n , n < , 



E/^latt 

k 



(ii/ n + \u m )G% xt (iv n )5 n!nl , (16) 



where 



G k a (iv n ) = [w n + n - e k - Ei oc (ii/ n )] 



(17) 



and Gioc, Sioc are the Green's function and self-energy 
of the local DMFT impurity problem, respectively. The 



equivalent susceptibilities in the mixed Legcndrc-Fouricr 
representation are simply obtained as the inverse of 
Eq. (14) 



X«'(iw m ,q)= 2J T niX°(ivn,iv n ',iuj m ,q)T n 



(18) 



where the high-frequency behavior of G\ oc (iu n ) and 
G latt (ii/ n ) can easily be considered in the frequency 
sums. Evaluation of lattice susceptibilities from 
Xiatt (i^n, hv j iw m ) can also directly be propagated to the 
mixed Legendrc-Fourier representation, abolishing alto- 
gether the need to transform back to Fourier representa- 
tion 



X(iw m ,q) = -g2 ^2 Xiatt(i^n,hvn',iw m ,q) 

& ^(-i) w 'V2TTTV27 7 TTxiattj/'(iw m ,q). 

(19) 



/3 2 



Note that xiatt can be written as the sum of a free two- 
particle propagation Xj att Eq. (16) (bubble part) and a 
connected rest xfat" 11 (vertex part). These two terms can 
be separately summed in Eq. (19). 

The present mixed basis representation has been suc- 
cessfully used in a recent investigation of static finite- 
temperature lattice charge and magnetic susceptibilities 
for the Na a: Co02 system at intcrmcdiatc-to-largcr dop- 
ing x. 2e A first example for the dynamical, i.e., finite- 
frequency, case will be discussed in Sec. II C. 



B. Antiferromagnetic susceptibility of the 
three-dimensional Hubbard model 



In order to benchmark our approach, we investigate the 
antiferromagnetic susceptibility of the half-filled Hub- 
bard model (3) on a cubic lattice within the DMFT 
framework. All quantities are again expressed in units 
of the hopping t and with U/t = 20 and T/t = 0.45. 
This temperature is sufficiently close to the DMFT Neel 
temperature Tn ~ 0.30t to yield a dominant vertex part, 
while still having a non-negligible bubble contribution. 

We compute the susceptibility xioc of the DMFT im- 
purity problem using the CT-HYB algorithm. In Fig. 7 
we compare the mixed Legendre-Fourier coefficients 
X;z'(iw m ) to the Fourier coefficients x(i^„, \v n ' , iuj m ). For 
clarity, we focus on the first bosonic frequency iui m = 0. 
We observe that the xw (0) have a very fast decay except 
in the I = V direction. This contrasts with the behavior 
of xi}v n , Wn' 1 0) which exhibits slower decay in the three 
major directions \v n = 0, \v n i = and \v n = \v n i . 

The generalized susceptibility in r-diffcrences (11) has 
discontinuities along the planes T14 = and T14 = t\i + 
T34 as well as non-analyticitics (kinks) for T12 = and 
T34 = 0. These planes induce corresponding slow decay 
in the Fourier representation (14). 22 When it comes to 
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FIG. 7. (Color online) Generalized local magnetic suscepti- 
bility x£c = |(XkI _ XiL) at tne bosonic frequency iu m = 
computed from the DMFT impurity problem. Upper panel: 
coefficients x™ c u , (0) in the mixed Legendre-Fourier represen- 
tation. Lower panel: Fourier coefficients xS c (ii / n, w n ' > 0)- 



the mixed Legendrc-Fouricr representation (13) however, 
the planes T\% = and r 34 = are on the border of 
the imaginary-time region being expanded in this basis, 
which renders the coefficients insensitive toward these. 

Computing lattice susceptibilities from Eq. (15), it 
is necessarily required to truncate the matrices. This 
leads to difficulties when computing the susceptibility 
from the Fourier coefficients xi oc (ii/„, i^„', icj m ). As we 
can see from Fig. 7, the Fourier coefficients have a 
slow decay along three directions. The inversion of 
Xioc(itVi, izVi iw m ) is delicate because many coefficients 
are involved even for large v 1 . One needs to use a very 
large cutoff to obtain a precise result. Alternatively, one 
can try to separate the high- and low-frequency parts of 
the equation and replace the susceptibilities with their 
asymptotic form at high frequency (see Ref. 27). While 
is is effectively possible to treat larger matrices, it is still 
required to impose a cutoff on the high-frequency part 
for the numerical computations. 

In the mixed Legendre-Fourier representation, the sit- 
uation is different. Only the coefficients along the di- 
agonal decay slowly. In the inversion of the matrix, the 
elements on the diagonal for large / are essentially recom- 
puted from themselves. One can expect that there will 
be a lot less mixing and thus a much faster convergence 
function of the truncation. 

In Fig. 8, we display the vertex part of the generalized 
lattice susceptibility Xiatt - Xiatt obtained from Eq. (15) in 
both representations. In both cases, we see that the diag- 
onal part quickly becomes very small. In other words, the 
diagonal of the lattice susceptibility is essentially given 
by the bubble part x?atr However, while essentially all 



0.40 r 




FIG. 8. (Color online) Vertex part of the generalized mag- 
netic lattice susceptibility x£tt — Xiatt a ^ the bosonic fre- 
quency iio m = and at the antiferromagnetic wave vector 
q = (7r,7r,7r). Upper panel: coefficients xStt.M' (°) -X?att,H' (°) 
in the mixed Legendre-Fourier representation. Lower panel: 
Fourier coefficients xSttO'V iu n >,0) - XiattO^n* «V,0) of the 
lattice susceptibility. Both plots employ the same number of 
coefficients. 
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FIG. 9. (Color online) Antiferromagnetic susceptibility as 
a function of the number of Legendre (#Z = / max + 1) and 
Matsubara (#n = 2n max + 2) coefficients, respectively, used 
in the calculation. 



the information is condensed close to 1,1' = in the 
mixed Legendrc-Fouricr representation, the Fourier coef- 
ficients still have a slow decay along the directions given 
by iv n = and \v n i = 0. From this figure one can specu- 
late that a quantity computed from the Legendre-Fourier 
coefficients will converge rapidly as a function of a cutoff 
l mM . However, we need to make sure that the coefficients 
close to I, I' = arc not affected much by the truncation. 

In order to assess the validity of these speculations 
we compute the static antiferromagnetic (q = (tt, 7T, 7t)) 
susceptibility x m (0, q) as a function of the cutoff in both 
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FIG. 10. (Color online) Imaginary part of the magnetic sus- 
ceptibility on the real frequency axis along high-symmetry 
lines in the 2D Brillouin zone. 



representations. It is obtained from Eq. (19) using the 
magnetic susceptibility x m = — x H ). 

Since the diagonal of the lattice susceptibility is essen- 
tially given by the bubble (see Fig. 8), the sums above are 
performed in two steps. The vertex part shown in Fig. 8 
is summed up to the chosen cutoff, while the bubble part 
is summed over all frequencies with the knowledge of its 
high-frequency behavior. The result is shown in Fig. 9. 
It reveals a major benefit of the Legendre representation: 
the susceptibility converges much faster as a function of 
the cutoff. The static susceptibility is essentially con- 
verged at l max ^12. This corroborates the idea that the 
small-Z, I' part of xiatt,H' is only weakly dependent on the 
further diagonal elements of Xioc.ii'- 



C. Dynamical susceptibility of the two-dimensional 
Hubbard model 

As a final benchmark, we demonstrate that our method 
is not restricted to the static case. To this end, we show 
the momentum resolved dynamical magnetic susceptibil- 
ity x(w,q) for a DMFT calculation for the half- filled 
two-dimensional (2D) square lattice Hubbard model in 
Fig. 10. We have chosen an on-site interaction Ujt = 4 
and temperature T/t = 0.25, which is slightly above 
the DMFT Neel temperature. The susceptibility was 
computed from the Legendre representation according 
to Eq. (19) using 20 x 20 Legendre coefficients, which 
was sufficient for all bosonic frequencies. In general, for 
higher bosonic frequencies more Legendre coefficients are 
needed to represent the vertex part of the generalized 
magnetic lattice susceptibility. However, no additional 
structure appears in the high I, I' region. We then ana- 
lytically continued the data using Pade approximants. 28 
The figure shows the typical magnon spectrum 6 ' 29,30 rem- 
iniscent of a spin wave in this paramagnetic state with 
strongly enhanced weight at the antiferromagnetic wave 
vector q = (n, it) due to the proximity of the mean-field 
antiferromagnetic instability. 



III. CONCLUSION 

In this paper, we have studied the representation of 
imaginary-time Green's functions in terms of a Legen- 
dre orthogonal polynomial basis. We have shown that 
CTQMC can directly accumulate the Green's function 
in this basis. This representation has several advan- 
tages over the standard Matsubara frequency represen- 
tation: (i) It is much more compact, i.e., coefficients de- 
cay much faster; this is particularly interesting for stor- 
ing and manipulating the two-particle Green's functions. 
Moreover, two-particle response functions can be com- 
puted directly in the Legendre representation, without 
the need to transform back to the Matsubara represen- 
tation. In particular, the matrix manipulations required 
for the solution of the Bethe-Salpeter equations can be 
performed in this basis. We have shown that this greatly 
enhances the accuracy of the calculations, since in con- 
trast to the Matsubara representation the error due to 
the truncation of the matrices becomes negligible, (ii) 
The Monte Carlo noise is mainly concentrated in the 
higher Legendre coefficients, the contribution of which 
is usually very small; this allows us to develop a system- 
atic method to filter out noise in physical quantities and 
to obtain more accurate values for, e.g., the correlation 
energy in LDA+DMFT computations. 
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Appendix A: SOME PROPERTIES OF THE 
LEGENDRE POLYNOMIALS 

In this appendix, we summarize for convenience some 
basic properties of the Legendre polynomials. Further 
references can be found in Rcfs. 31, 33, and 34. We use 
the standardized polynomial Pi(x) defined on x £ [—1, 1] 
through the recursive relation 

(I + l)P l+ i(x) = (21 + l)xl\(x) - lPi-i{x) (Al) 
P (x)=l, P 1 {x)=x (A2) 
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Pi arc orthogonal and their normalization is given by 



dxP k (x)Pi (x) = T^—Ski (A3) 



The Pi are bounded on the segment [—1, 1] by 34 

\Pi(x)\ < 1 (A4) 

with the special points 

Pi(±i) = (±iy. (as) 

The primitive of Pi(x) that vanishes at x = —1 is (cf. 
Ref. 31, Vol II, section 10.10) 



J &yPi(y) = 



P+i{x) - Pl-^X) 
21 + 1 



I > 1 



(A6) 



By orthogonality or (A5), it also vanishes at x = 1. 
The Fourier transform of the Legendre polynomial re- 
stricted to the segment [—1, 1] is given by formula 7.243.5 
of Ref. 33 

f^Pi(x)dx = i l ^J l+h {a) 

= 2i l j l {a), (A7) 

where J denotes the Bessel function and ji(a) — 
\p52<h+\ (°) denotes the spherical Bessel functions. 



Appendix B: FAST DECAY OF THE LEGENDRE 
COEFFICIENTS 



Let us consider a function g[r) smooth on the seg- 
ment [0,/3] (i.e. to be precise C°°, indefinitely diffcrcn- 
tiable), and j3— antiperiodic, like a Green's function. In 
this appendix, we show that its Legendre coefficients de- 
cay faster than any power law contrary to its standard 
Fourier expansion coefficients which decay as power laws 
determined by the discontinuities of the function and its 
derivatives. 

Let us start by reminding the asymptotics of the stan- 
dard Fourier expansion coefficients on fermionic Matsu- 
bara frequencies. These coefficients are given by 



g(iv n )= f drg{T)e^ T 
Jo 



dr g'{r) 



(Bl) 



(B2) 



The coefficients vanish for n — > oo , and applying the same 
result to g' , one obtains 



s .g(/n + ff(0 + ) , n f 1 
g{iv n ) = : h O 



(B3) 



a function f(x) smooth on [—1, 1]. We can proceed in a 
similar way using the primitive of the Legendre polyno- 
mial (which is also given by a simple formula, (A6)). For 
I > 1, we have 



fl 



V21 + T 



dx f{x)Pi(x) 



/(*) 



dy Pi{y) 



^dx f{x) (J dyPi{y) 



The crucial difference with the Fourier case is that, for 
I > 1, the boundary terms always cancel, whatever the 
function / due to the orthogonality property of the poly- 
nomials (it can also be checked directly from (A5)). So 
we are left with just the integral term. Since the Legen- 
dre coefficients of /' vanish at large I (by applying the 
previous formula to /'), we get instead of (B3) 



fl 



{ I 



(B5) 



In both cases, the reasoning can be reproduced recur- 
sively, by further differentiating the function, as long 
as no singularity are encountered. In the Fourier case, 
it produces the well-known high-frequency expansion in 
terms of the discontinuity of the function and its deriva- 
tives. In the Legendre case, we find that the coefficients 
arc o{l/l k ) as soon as / is k times differentiable. Hence if 
the function is smooth on [—1, 1], the coefficients decays 
asymptotically faster than any power law. 

The only point that remains to be checked is that in- 
deed G(t) is smooth on [0, /?]. It is clear from its spectral 
representation 



G{t) 



dv- 



1 



—AM 



(B6) 



if we admit that the spectral function A{v) has compact 
support, by differentiating under the integral. 

Finally, while this simple result of "fast decay" is 
enough for our purposes in this paper, it is possible to 
get much more refined statements on the asymptotics of 
the Legendre coefficients of the function /, in particular 
when it has some analyticity properties. For a detailed 
discussion of these issues, and in particular of the condi- 
tions needed to get the generic exponential decay of the 
coefficients, we refer to Ref. 35. 



Appendix C: DIRECT ACCUMULATION OF THE 
LEGENDRE COEFFICIENTS FOR THE CT-HYB 
ALGORITHM 



Let us now turn to the Legendre expansion. Using the 
same rescaling as before, we can consider for simplicity 



In this appendix, we describe how to compute directly 
the Legendre expansion of the one-particle and the two- 
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particle Green's function. 



1. The accumulation formulas in CT-HYB 

For completeness, let us first recall the accumulation 
formula for the one-particle and the two-particle Green's 
functions in the CT-HYB algorithm 8,9,11,19 , which sums 
the perturbation theory in the hybridization function 
A ab (\v n ) on the Matsubara axis. While these formulas 
have appeared previously in the literature, this simple 
functional derivation emphasizes the "Wick" -like form of 
the high-order correlation function. 

The partition function of the impurity model reads 



Z = J Pc t 2?ccxp(-S' cff ) 



(CI) 



where the effective action has the form 



Sett 



-dr'J2c\(r)G-^ B (T,T')c B (T') 

A,B 



+ / dTH int ({ C \(T),C A (T)}) 



(C2) 



G 0AB( ilJ n) = (i^™ + I^&ab - h° AB - A AB (iv n ), (C3) 

To simplify the notations, we use here a generic index 
A, B. In the case where there are symmetries, like the 
spin SU(2) symmetry in the standard DMFT problem, 
the Green's functions are block diagonal. For example, 
the generic index A can be (a, a), where a is an orbital 
or site index, and spin index a ^ is the block index. 

The partition function is expanded in powers of the 
hybridization A as 

/n 
J] dndri £ w(n, {X 3 , \' v Tj , rj}) (C4) 

nxu i=l A;, A; 

w(n, {Xj, Aj, Tj, Tj}) = ^e^jA^y.in - rj)] x 

Trfre-^flcliT^Tl)), (C5) 



(=i 



where T is time ordering and H\ oc is the local 
Hamiltonian 8,9,11,19 . \w\ are the weights of the Quan- 
tum Monte Carlo Markov chain. Introducing the short 
notation C = (n,{Xj,X'j,Tj,Tj}) for the QMC configu- 
ration, the partition function Z and the average of any 
function / over the configuration space (denoted by an- 
gular bracket in this section) are given by 



Z = J>(C) 
c 

(/(C)) = i£>(C)/(C) 



(C6) 
(C7) 



The one-particle and two-particle Green's functions are 
obtained as functional derivatives of Z with respect to 
the hybridization function, as 



Gab(ti,t 2 ) 



G ABCD (Tl,T2,T 3 ,T i ) 



1 dZ 
'ZdAgAfan) 
1 d 2 Z 



(C8a) 



Z d A B a , n )d A DC (r 4 , r 3 ) 
(C8b) 



In order to use the expansion of Z, we need to com- 
pute the derivative of a determinant with respect to its 
elements. Let us consider a general matrix A, its inverse 
M = A -1 and use Grassman integral representation 



dct A 



'[drjidrji 



(C9) 



Using the Wick theorem, we have 
d dct A 



dA 



ba 



~\dr}idfji(f\ b r) a )e^ * Ay% 



det A x M ab 



(ClOa) 



d 2 dct A 
dA ba dA dc 



det A (M ab M cd - M ad M cb ) (ClOb) 



Let us now apply (CIO) by introducing for each configu- 
ration C = (n, {Xj, X'j, Tj, t^}) the matrix A(C) of size n 
given by 

AiQij^A^y.in-T'j) (Cii) 
and its inverse M c = (A(C)) . We obtain 



dw(C) 



'(C) 



ddctA(C) dA(C)p a 
dA BA (T 2 ,Ti) detA(C) dA{C)f} a OAba^tx) 



E 



(r , M c dA{C) Pa 

W {C) 2^ M a pjTT 7— — T 

aT-l P OA BA {T 2 ,Tl) 



(C12a) 



and 



d 2 w{C) 



;(C) 



dA BA (T 2 ,Ti)dA DC {T4,,T S ) detA(C) 
y> d 2 dctA(C) dA(C)p a dA{C) Sl 
J^ =1 dA{C) Pa dA(C) 8l dA BA (T 2 ,n) dA DC ( Ti ,T 3 ) 



(C12b) 



Denoting 



D(C)f BTiT2 



dA(C) 0c 



dA BA ( T2 ,n) 

<>{ T 1 _ T 'a)H T 2 ~ Tp)6\' at A6\ p ,B (C13) 
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we finally obtain the accumulation formulas for the 
Green's functions 8,9 



Finally, Eq. (C18) can be transformed to a measurement 
in the Legendre representation according to (2) 



G AB (n,r 2 ) = - ( E M apD(C)f BTiT2 ) (C14a) 

W=i / 

/ ™ 

Gi 4 L^i,r 2 ,T 3 ,T 4 ) = / £ {M c ap M c lS -M c a5 M c lfi )x 



D(cr/ BT1T2 D(cr c ° DT3 



7<5 



(C14b) 



2. Legendre expansion of the one particle Green's 
function 



We take into account the time translation invariance 
and the r-antiperiodicity of the Green's function in the 
following way. A priori, in (C14a), the arguments ti,t 2 
are in the interval [0,(3]. We can however easily make 
this function j3— antiperiodic in both arguments 

Gab(ti, t 2 ) = 

- ( £ M ^~(n - r' a )5-{r 2 - T P )5 K ^ A 5 Xfj , B \ 

(C15) 

where we defined the periodic and antiperiodic Dirac 
comb respectively by 



= £(±l)»<f(r-n0) 



(C16) 



At convergence of the Monte-Carlo Markov chain, the 
Green's function is in fact translationally invariant in 
imaginary time and we have 



1 r 

G A b(t) = - dsG AB (T + s,s) 
P Jo 



(C17) 



which leads to 



/ ™ 

Gab{t) = — ( J2 M apt-{r « Tp))8 K , A 8 Xp , 

P W=l 



7 

(C18) 



P \a/3=l 

where P(St) is defined by 



Pi(8t) 



Pi(x(St)) St>0 
-Pi{x{5t + P)) 5t<0 



(C19) 



(C20) 



3. Legendre accumulation of the two-particle 
Green's function 



The generalized susceptibility \ °f (H) can be ex- 
pressed in term of G and G^ as 



Xabcd( T 12 , T34, Tl 4 ) =G { ^ a(y d(y , C(y , (t 2 1 , T 43 , T 23 ) 

— Gb tT)QCT (T 2 i)G £ i CT ' iC(T '(T43) (C21) 



so in this subsection we will focus on the computa- 
tion of G' 4 -*. We take into account the time transla- 
tion invariance with the same technique as for the one- 
particle Green's function. First we make the function 
G( 4 '(ti, r 2 , r 3 , r 4 ) fully (3— antiperiodic in the four vari- 
ables using the antiperiodic Dirac comb S~ defined in 
(C16), and we use the time translation invariance of the 
Green's function to obtain 



G (4) (ti 2 ,t 34 ,ti 4 ) 
If 13 

= — I df G (4) ("T14 + f , T14 - T 42 + f , T 34 + f , f) 

P Jo 

(C22) 

From (C14b), we get 



J 



Gf BCD ^ , T34, r 14 ) = -l £ (MZpMfs ~ MZ 5 M^)X 



5-(t 12 - (r' a - rp)) 8-(t 34 - (r; - r s )) 6+ (r 14 - (r' a - r 5 )) <^,a<^,b <^A' ,C^A 5) d ) (C23) 
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where S + and 6 arc defined in (C16). 
is straightforwardly obtained as 



Applying (13), the accumulation formula in the mixed Legendrc-Fouricr basis 



V2T+1V2F + 1 



(-1) 



i'+l. 



/ ( M Zp M 's ~ MSsM%p)Pl {r' a - T p) Pv « - n) e lu ^<- T ^8 K ,A5x„B8 K ,c5x s ,D ) (C24) 



where P is defined in (C20). 



We note that the measurement can be factorized to 
speed up the measurement process. In the Legendre mea- 
surement, only the part involving the first product of M- 
matrices factorizes, as can be seen from (C24). Note, 
however, that the second product of M-matrices merely 
generates crossing symmetry, so that the full informa- 
tion on this quantity is already contained in the first 
term. Hence this symmetry can be reconstructed after 
the simulation. In the one band case, the second product 
is proportional to 5 aa i , so that the G^^-component can 
be measured directly. For the G^ 4 ^^-component we only 
measure the term proportional to M^M^ and construct 
this component by antisymmetrization afterwards. 



Appendix D: ACCUMULATION FORMULA FOR 
THE CT-INT AND CT-AUX ALGORITHMS 



10°^ 
10-' 

io- 2 
io- 3 

E?10- 4 

io- 5 

1(T 6 

io- 7 
io- 8 



FIG. 11. (Color online) \T n i \ for the first even (red) and odd 
(blue) Legendre coefficients. The high-frequency tail is repro- 
duced correctly by t\ p \ 




Using a notation in analogy to the previous section, 
the expansion of the partition function Z, Eqs. (C1-C3), 
in the continuous-time interaction expansion (CT-INT) 
method 7 is given by 

/n 
Y[dT t £ wfaiXj^Tj}) (Dl) 

n^u t=l X 2i -i,X' 2i _ 1 
X 2 i,X 2i 

w(n, {Xj , Xj, Tj\) = —7 det [G 0Ai y (n - fj )] x 

n 

1=1 

where fj = T|(»+i)/2J an( i we have assumed the in- 
teraction part of the Hamiltonian to be of the form 
Hint({c\,c A }) = Y,abcd u abcdc\c b c { c c d and A = 
(a, a) is a generic index with a being the orbital or site 
index and a =t,| the spin index. In the CT-INT algo- 
rithm, we propose to measure the Legendre coefficients of 
S = SG based on the self-energy binning measurement 
originally introduced for the continuous-time auxiliary 
field (CT-AUX) algorithm 10 . Introducing the matrix 

G (C)„- =G 0Ai v.(^-^) (D3) 



and its inverse, M c = (Gq(C)) l , the self-energy binning 
measurement for the CT-INT can be written as 



Sab(t) = - ( ]T S(t - f a )S AK M^G°^ B (f ) \ . 

W=i / 

(D4) 

This can be straightforwardly transformed to a measure- 
ment in the Legendre basis by applying (2): 



S A b,i = -v / 27Tl( ]T ^P(*(r a ))M^G A3B fo)). 

\a/3=l I 

(D5) 

An analogous formula also applies to the CT-AUX. In 
practice, translational invariance may be used to gener- 
ate multiple estimates for S within a given configuration. 
The Green's function is obtained by transforming S to 
Matsubara representation and using Dyson's equation. 
The moments of G are straightforwardly computed from 
the moments of EG and the knowledge of those of Gq. 
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Appendix E: EXPLICIT FORMULA FOR T n! AND 
ITS HIGH FREQUENCY EXPANSION 



We will now deduce the coefficients of the expan- 
sion of T n i 



The transformation matrix from the Legendre to the 
Matsubara representation is 

Tm = ^+1 I dre^P (x(t)) (El) 
P Jo 

where v n is a fermionic Matsubara frequency and I the 
Legendre index. Using (A7) and introducing the reduced 
frequencies D n = fiv n = (2n + 1)7T, we find 



Tm = (-1)" i ,+1 V2r+lj« (y) 



(E2) 



Note that T n i is actually independent of /?. 

T n i is a unitary transformation, as can be check ex- 
plicitely using the Poisson summation formula and the 
orthogonality of the Legendre Polynomials (A3) 



An) 



P >l y nJ 



(E4) 



^ V2TTWWTT 



f f drdr' PMr))Pv {x(r')) - ]T e -^-r') 
J Jo , 



=5(t-t') 

--V2TTTV2TTT / — Pi{x)p v {x) 

-S w (E3) 



This straightforwardly done from an corresponding repre- 
sentation of the Bcssel function, cf. e.g. Rcf. 36, Section 
10.1 



Li 



ji(z) = z- 1 x ^ sin(z - 7rf/2) 



lk (l + W(2z 



,-2k 



k=0 



{2k)\{l - 2k)\ 



-cos(z-7r//2) ^ (-1) 



k=0 



fc (/ + 2fc+l)!(2z)- 2fc - 1 
(2k + l)\(l - 2k- 1)! 



(E5) 



For the case at hand this gives 



i 2v / 2TTT { cos )Y J 



Li 



(l + 2k)\ 



k=0 



{2k)\{l-2k)\ {W n ) 



2k+l 



(l + 2k + l)\ 



1 



fe=0 



(2* + — 2fc — 1)! (iD n ) 2k + 2 



(E6) 



r 



The two sums can be combined to 



Tnl = 2V2TTTJ27- A j 

p—1 ^ 



(i+ P -iy. (-iy 



=J(P-I)!(l-P+I)!(ii?n)" 



^p+;,odd, 



(E7) 

which immediately provides the coefficients of (E4) 



(-l) p 2 v / 2T+T 



(i+p-iy. 

(p-iy.Q-p + iy 



Jp+l, odd- 



(E8) 



Fig. 11 shows T n i for the first Legendre coefficients plot- 
ted against the fermionic Matsubara frequency \v n . The 
doubly logarithmic plot clearly shows the high-frequency 
l/i^„-behavior for the even and the l/(ii/ n ) 2 -behavior for 
the odd coefficients. One can see that, as expected, struc- 
ture at very high frequencies is only carried by polyno- 
mials with large values of I. 
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